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Abstract 


The  use  of  finite-difference  and  finite-element  computer  codes  to  solve  problems  involving 
fast,  transient  loading  is  commonplace.  A  large  number  of  commercial  codes  exist  and  are 
applied  to  problems  ranging  from  fairly  low  to  extremely  high  damage  levels  (e.g.,  design  of 
containment  structures  to  mitigate  effects  of  industrial  accidents;  protection  of  buildings  and 
people  from  blast  and  impact  loading;  foreign-object  impact  damage;  and  design  of  space 
structures  to  withstand  impacts  of  small  particles  moving  at  hypervelocity,  a  case  where  the 
pressures  generated  exceed  material  strength  by  an  order  of  magnitude).  But,  what  happens  if 
code  predictions  do  not  correspond  with  reality?  This  report  discusses  various  factors  related 
to  the  computational  mesh  that  can  lead  to  disagreement  between  computations  and  experience. 
Subsequent  reports  will  focus  on  problems  associated  with  contact  surfaces  and  material 
transport  algorithms,  constitutive  models,  and  the  use  of  material  data  at  strain  rates 
inappropriate  to  the  problem.  It  is  limited  to  problems  involving  fast,  transient  loading,  which 
can  be  addressed  by  commercial  finite-difference  and  finite-element  computer  codes. 

This  report  has  been  accepted  for  publication  in  a  future  volume  of  the  International  Journal 
of  Impact  Engineering. 
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1.  Introduction 


This  report  focuses  on  the  numerical  simulation  of  problems  in  mechanics  involving  fast, 
transient  loading.  For  practical  purposes,  these  can  be  divided  into  problems  involving  structural 
dynamics  and  wave  propagation  (Table  1).  There  is  no  clear  demarcation  between  these  two  areas. 
The  labels  are  misleading  since  both  types  of  problems  involve  wave  propagation.  Nonetheless, 
these  designations  have  caught  on  in  the  literature.  These  two  labels  deal  with  the  behavior  of  inert 
materials  that  are  subjected  to  intense  impulsive  (distributed  over  a  surface,  such  as  air  blast,  over 
a  long  time  [milliseconds  to  seconds])  or  impact  (applied  to  a  single  point  or  a  very  small  area  over 
a  very  short  [nanoseconds  to  microseconds]  time  span)  loading. 

There  is  also  a  large  class  of  energetic  materials  that  reacts  quite  differently  when  excited. 
Energetic  materials  are  not  discussed  in  this  report.  For  more  information  on  this  topic  see  Zukas 
and  Walters  [1],  Cooper  [2],  Cooper  and  Kurowski  [3],  Cheret  [4],  Blazynski  [5],  Fickett  and 
Davis  [6],  and  Mader  [7, 8]. 

Most  of  the  work  done  in  the  area  of  fast,  transient  loading  is  experimental  in  nature.  This  is 
due  either  to  complexities  of  geometry  or  the  nonlinearity  of  material  behavior  or  both.  Closed-form 
analytical  solutions  are  generally  rare  and  apply  only  to  some  small  subset  of  the  overall  problem. 

Numerical  solutions,  in  the  form  of  finite-difference  and  finite-element  codes,  have  been 
successfully  used  in  the  past.  In  particular,  the  combination  of  experiments,  numerical  solutions,  and 
dynamic  material  characterization  has  been  shown  to  be  very  effective  in  reducing  both  manpower 
requirements  and  cost  [9].  However,  the  computer  codes  available  for  dynamic  analyses  are  quite 
complex.  Considerable  experience  with  both  the  codes  and  the  physical  problems  they  are  intended 
to  solve  is  vital.  Also  critical  is  the  determination  of  material  constants  for  the  various  constitutive 
models  available  at  strain  rates  appropriate  to  the  problem. 
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Table  1.  Dynamic  Situations 


STRUCTURAL  DYNAMICS 

•  Blast/Shock  Loading  of  Structures 

•  Underwater  Explosions 

•  Fluid-Structure  Interactions 

•  Mechanical  System  Dynamics 

-  machinery  &  mechanisms 

-  agricultural,  construction,  off-highway  equipment 

-  turbomachinery  systems 

-  containment  structures 

-  vehicular  collisions 

-  aeronautical/aerospace  systems 

•  Biodynamic  Systems 

•  Plate  and  Shell  Structures 

•  Nonperforating  Impacts 

•  Rotating  Machinery 

•  Metal  Forming 

WAVE  PROPAGATION 

•  Lunar/Planetary  Impact 

•  Explosive  Welding,  Forming,  Compaction 

•  Shock  Consolidation/Shock  Synthesis 

•  Chemical  Energy  Penetrators 

-  explosively  formed  penetrators  (EFP) 

-  shaped  charge  jets 

•  Kinetic  Energy  Penetration 

-  fragments 

-  long  rods 

-  bombs 
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Figure  1  shows  the  dramatic  effect  that  experience  can  have  on  computed  results.  The  figure, 
provided  by  Dr.  Paul  Senseny  of  the  former  Defense  Special  Weapons  Agency  (now  called  Defense 
Threat  Reduction  Agency),  shows  the  results  obtained  by  four  different  users  of  the  DYNA  code  for 
a  problem  involving  airblast  loading  of  a  silo  door.  Each  of  the  code  users  worked  independently. 
Note  that  four  people  with  four  different  backgrounds  produced  four  distinctly  different  results  with 
the  same  code  in  solving  the  same  problem. 

Why  this  emphasis  on  experience?  To  be  sure,  experience  is  not  a  sufficient  factor  to  guarantee 
accurate  computations.  Indeed,  all  four  code  users  in  the  previous  problem  could  have  been 
experienced  in  some  sense.  But  experience  is  an  important  factor.  There  exist  many  cases  in 
computational  continuum  dynamics  where  tradeoffs  must  be  made  to  achieve  reasonable  results  in 
a  reasonable  time  at  finite  costs.  Szabo  and  Actis  [10]  point  out  the  importance  of  timely  solutions 
in  industry.  Quoting  from  an  engineer’s  experience  at  a  major  industrial  laboratory,  they  point  out 
that  “...the  stress  engineer  may  not  be  able  to  guide  the  design  because  by  the  time  he  has  generated 
his  several  thousand  degree-of-freedom  finite-element  model,  prototypes  are  already  being  made.” 
Anyone  can  perform  finely  resolved  one-dimensional  (1-D)  calculations  almost  by  rote. 
Two-dimensional  calculations  (2-D),  though  now  largely  performed  on  workstations  and  personal 
computers,  require  a  keen  knowledge  of  the  problem  and  numerical  simulation  methods.  This  holds 
even  more  for  successful  three-dimensional  (3-D)  simulations  that,  with  all  their  compromises,  will 
still  be  expensive  and  require  significant  computer  resources. 

2.  Difficulties  Inherent  in  Numerical  Approximations 

When  a  user  acquires  a  code,  he/she  receives  a  package  in  which  certain  decisions  have  been 
made  by  the  developer.  The  user  has  no  control  over  some  of  the  approximations  that  come  as  part 
of  the  package  but  must  be  aware  of  them,  as  they  will  affect  all  solutions  that  the  user  will  generate 
with  the  package.  This  section  briefly  describes  these  errors  and  cites  sources  in  the 
finite-difference/finite-element  literature  where  extensive  discussions  can  be  found.  These  are 
problems  inherent  in  transforming  a  physical  problem  into  a  discrete  model  and  solving  it  on 
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Displacement,  cm 


The  Same  Code  Will  Produce  Different  Answers 
In  The  Hands  of  Different  Users 


Time,  psec 


DYNA  was  used  in  all  calculations 


Figure  1.  Variation  of  Computational  Results  With  User  Experience. 


computers  with  limited  precision  through  the  medium  of  finite  differences  and  finite  elements.  Each 
difference  scheme  and  each  element  have  unique  characteristics  that  can  affect  the  numerical 
solution.  For  example,  most  implicit  schemes  for  advancing  the  solution  in  time  have  a  certain 
amount  of  viscosity  built  into  them  over  which  the  user  generally  has  no  control.  Yet,  the  user  must 
be  aware  of  this  characteristic  when  he/she  sees  his/her  numerical  solution  drift  out  of  phase  with 
an  experimental  or  analytical  result.  Certain  finite  elements  are  prone  to  locking  (becoming 
excessively  stiff)  when  the  constitutive  relationship  is  evaluated  at  element  Gauss  points.  This  can 
be  overcome  by  evaluating  the  constitutive  equations  at  one  point  (usually  the  element  centroid). 
This  procedure  is  called  “underintegration.”  It  eliminates  locking  but  gives  rise  to  spurious 
deformation  modes  known  as  “hourglassing.”  A  user,  unaware  of  this,  might  interpret  such 
numerical  noise  as  a  physical  response.  Explicit  integration  schemes  are  only  conditionally  stable. 
Using  stability  fractions  built  into  explicit  hydrocodes  without  consideration  of  the  problem  being 
solved  can  lead  to  unstable  solutions  that  can  degrade  over  hundreds  of  cycles  and  could  be 
interpreted  by  novice  code  users  as  a  physical  response  of  the  system  being  modeled. 

The  errors  inherent  in  finite-element  modeling  are  nicely  discussed  by  Utku  and  Melosh  [11], 
as  are  guidelines  for  mesh  preparation  by  Melosh  and  Utku  [  1 2] .  A  superb  work  on  practical  aspects 
for  static  and  dynamic  finite-element  analyses  is  the  book  by  Meyer  [13].  Extremely  valuable 
insights  into  finite-element  design  and  its  practical  application  are  to  be  found  in  MacNeal’s 
book  [14].  Morris  and  Vignjevic  [15]  review  error  control  and  error  bounding  methods  for  finite 
elements.  They  also  present  a  method  for  error  control  in  the  idealization  phase  of  a  full 
finite-element  analysis.  These,  together  with  personal  experience,  are  the  main  references  for  the 
material  in  this  section. 

An  analyst  begins  by  examining  some  physical  problem  that  needs  to  be  modeled  numerically. 
This  may  be  necessary  because  it  is  too  expensive,  too  difficult,  or  just  plain  impossible  (as  in 
nuclear  testing)  to  perform  parametric  experiments.  It  might  be  because  the  information  needed  does 
not  lend  itself  to  direct  measurement.  It  is  often  the  case  [9,  16]  that  the  combined  use  of 
experiments  and  calculations,  including  some  material  characterization  at  high  strain  rates,  produces 
more  information  in  less  time  and  at  less  cost  than  reliance  on  experiments  or  calculations  alone. 
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At  first  glance,  a  problem  may  appear  to  be  extremely  complex.  Pressure  vessels  come  with  various 
cutouts,  end  caps  that  can  range  from  circular  plates  through  hemispherical  shells,  pipes  coming  in 
and  out,  some  irregular  sections,  weld  points,  and  a  countable  infinity  of  bolts  and  assorted 
restraining  devices.  Face  it,  real  structures  are  messy  to  model  numerically.  They  can  be  mounted 
to  floors  and  walls  and  receive  and  transmit  loads  to  these  through  a  variety  of  mounts  (shock 
isolation  is  popular  in  Japan  and  other  earthquake-prone  areas).  The  space  shuttle  has  a  nice,  clean 
shape,  but  most  space  structures  tend  to  show  greater  similarity  to  the  highly  irregular  surfaces  of 
the  space  cruisers  in  Star  Wars  than  to  the  shuttle.  Thus,  high-velocity  impact-generated  debris  and 
ricocheting  projectiles  can  interact  with  critical,  externally  mounted  equipment  such  as  antennas  and 
solar  panels. 

Three-dimensional  impact  problems  are  the  rule  rather  than  the  exception.  Uniform  grid 
resolution  is  generally  not  possible  in  practical  problems.  Some  gradation  of  the  mesh  is  required. 
If  this  is  not  done  with  great  care,  spurious  signals  and  assorted  numerical  artifacts  arise  in  the 
calculation,  leading  to  instabilities  or  a  masking  of  the  desired  results.  Even  coarse  calculations  can 
require  in  excess  of  40  central  processing  unit  (CPU)  hours  on  modem  computers.  If  sufficient 
memory  is  not  available  to  run  the  problem  “in-core,”  then  extensive  buffering  between  main 
memory  and  mass  storage  is  required,  further  increasing  turnaround  time  and  cost.  For  sufficiently 
large  problems  on  a  small  central  memory  machine,  it  is  possible  to  approach  situations  where  the 
bulk  of  CPU  time  is  spent  on  input-output  operations  and  only  a  small  fraction  is  spent  in  advancing 
the  solution,  rendering  the  computation  uneconomical. 

Does  this  mean  that  3-D  calculations  cannot  be  done  today?  Not  at  all.  Quite  a  few  practical 
problems  have  been  successfully  addressed  with  3-D  codes.  However,  compromises  are  required, 
and  these,  in  turn,  require  a  keen  understanding  of  the  physical  problem,  the  effects  of  discretization 
on  that  problem,  and  the  effects  that  numerical  artifacts  (such  as  uneven  resolution  in  different 
coordinate  directions,  mixing  of  implicit  and  explicit  integration  schemes  or  explicit-explicit 
partitions,  choice  of  mesh  or  element  type,  effects  of  sliding  surfaces  or  interfaces,  and  the  use  of 
various  viscosities  to  stabilize  computations)  have  on  the  solution. 
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All  the  details  and  materials  in  a  real  structure  usually  cannot  be  accounted  for  in  a  numerical 
simulation.  Hence,  the  analyst  must  now  make  an  analytically  tractable  model  without  sacrificing 
the  essential  elements  that  make  up  the  response  of  the  physical  structure.  The  first  order  of  business 
is  to  simplify  the  physical  system  by  taking  all  essential  geometric  and  material  features  that  govern 
its  response  into  account.  Some  simplification  of  geometry  and  lumping  of  masses  is  inevitable. 
Depending  on  the  situation,  one  may  be  able  to  employ  specialized  mathematical  models  such  as 
beam,  plate,  or  shell  theory  (or  a  combination  of  these).  Preliminary  analyses  need  to  be  made  to 
determine  whether  large  strains  and  rotations  constitute  a  part  of  the  response  or  whether  a  linear 
analysis  will  suffice.  A  number  of  uncertainties  are  introduced  in  this  process,  some  of  which  cannot 
be  quantified.  Other  uncertainties  are  due  to  variabilities  in  material  properties,  loading,  fabrication, 
and  other  factors.  For  example,  rolled  homogeneous  armor  (RHA)  steel  is  used  extensively  in 
military  construction.  It  can  safely  be  said  that  RHA  is  rolled.  It  is  also  used  as  armor.  However, 
the  military  specifications  that  govern  the  production  of  RHA  have  wide  tolerances  so  that  it  is 
anything  but  homogeneous.  Material  properties  (primarily  hardness)  are  known  to  vary  by  as  much 
as  10%  within  a  lot  of  RHA  and  up  to  30%  from  lot  to  lot.  This  makes  single  tests  (the  famous 
“one-shot  statistics”)  useless  and  correlation  between  numerical  results  and  experiments  unlikely 
unless  a  statistically  meaningful  number  of  tests  have  been  done.  Simple  go/no-go  ballistic  tests  can 
cost  upward  of  $2,000  each.  Instrumented  field  tests  can  run  from  $10,000  to  $100,000  each.  As 
a  rule,  then,  a  statistically  significant  data  set  is  almost  never  available. 

3.  Idealization 

The  ultimate  goal  of  this  idealization  process — the  transition  from  a  complex  physical  model 
to  a  simpler  one  but  incorporating  all  the  relevant  physics — is  a  mathematical  model  consisting  of 
a  number  of  equations  that  closely  represent  the  behavior  of  the  physical  model.  A  formal  seven-step 
process  was  proposed  by  Morris  and  Vignjevic  [15].  An  experienced  analyst  will  go  through  the 
procedure  guided  by  a  few  principles  and  much  insight  garnered  over  a  long  career.  The  time 
required  might  take  weeks  to  months,  depending  on  the  complexity  of  the  physical  system  and  the 
accuracy  required  in  the  analysis.  If  the  assumptions  of  the  mathematical  model  are  reasonable,  very 
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little  mathematical  modeling  error  results.  If  this  is  poorly  done,  say  if  the  height- to- span  ratio  of 
beam  theory  is  violated,  the  thickness-to-radius  ratio  for  thin  shells  is  not  satisfied,  a  poor  choice  of 
material  properties  is  made,  or  the  constitutive  model  omits  a  critical  item  such  as  dynamic  failure, 
serious  errors  can  be  incurred.  Such  synthesis  is  not  taught  in  schools  but  learned  in  apprenticeship 
with  an  experienced  modeler.  Fortunately,  in  the  vast  majority  of  cases,  this  is  done  very  well  so  that 
mathematical  modeling  errors  are  negligibly  small,  or  at  least  smaller  than  errors  committed  by  code 
users,  the  topic  of  the  next  section. 

A  solution  is  needed  for  the  mathematical  model.  Since  most  problems  involving 
high-strain-rate  loading  are  not  analytically  tractable,  recourse  is  made  to  a  computer.  In  the  process 
of  computing,  especially  in  matrix  operations,  situations  occur  where  differences  between  numbers 
of  almost  equal  magnitude  must  be  taken.  This  can  lead  to  situations  where  the  roundoff  error  can 
completely  overwhelm  the  computed  quantity.  Given  a  machine,  there  is  always  a  limit  to  the  mesh 
refinement  beyond  which  computed  quantities  may  be  100%  erroneous  [11].  Roundoff  errors  may 
be  kept  negligibly  small  by  using  “...longer  wordlength  machines  and  double  precision  arithmetic 
with  not  too  refined  finite  element  meshes”  [11]. 

Most  dynamic  analyses  for  problems  involving  wave  propagation  are  done  with  the  simplest 
elements — constant  or  linear  strain  triangles  and  quadrilaterals  in  2-D  computations,  tetrahedra  and 
hexahedra  in  3-D  computations.  Many  early  (mid-70s)  finite-element  codes,  such  as  Lawrence 
Livermore  National  Laboratory  ’  s  DYNA  and  DYSMAS-L,  incorporated  high-order  elements  in  their 
initial  formulations.  With  experience,  these  were  dropped.  Many  transient  response  calculations 
in  the  wave  propagation  regime  involve  the  presence  of  steep  stress  gradients  and  shock  waves.  It 
has  been  found  with  experience  that  there  is  marginal  increase  in  accuracy  but  considerable  increase 
in  cost  in  trying  to  model  what  are  essentially  discontinuities  with  higher  order  polynomials.  The 
characteristics  of  these  elements,  and  the  problems  of  locking  and  reduced  integration  associated 
with  them,  are  lucidly  discussed  in  the  book  by  MacNeal  [14].  The  original  literature  and  the  code 
manuals  should  also  be  reviewed  in  order  that  these  elements  not  be  misused. 
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Code  users  should  also  be  aware  of  the  viscosity  built  into  implicit  temporal  integrators  and  the 
conditional  stability  of  explicit  integrators.  Additional  information  can  be  found  in  Belytschko  and 
Hughes  [17],  Donea  [18],  and  Zienkiewicz  [19]. 

4.  The  Human  Factor 

It  is  generally  accepted  that,  if  the  idealization  from  physical  system  to  mathematical  model  is 
done  well,  the  errors  due  to  truncation,  roundoff,  and  other  properties  of  finite-element  or 
finite-difference  schemes  can  be  readily  detected  and  contribute  to  no  more  than  about  5%  of  the 
total  solution  error.  Solutions,  however,  can  be  totally  invalidated  by  poor  choice  of  mesh,  failure 
to  include  relevant  physics  in  the  constitutive  description,  using  a  limited  or  inappropriate  database 
to  evaluate  the  constants  of  a  constitutive  model,  failure  to  recognize  instabilities,  or  the  effects  of 
contact  surfaces  on  numerical  solutions.  In  short,  computational  techniques  come  with  certain 
built-in  limitations  that  are  easily  recognizable  and,  with  few  exceptions,  controllable.  To  really 
mess  up  requires  a  human. 

Some,  but  hardly  all,  of  the  errors  in  the  application  of  computer  codes  to  practical  problems 
involve  the  following. 

Meshing:  A  code  user  has  available  a  wide  choice  of  elements  in  any  commercial  code. 
Having  selected  one  or  more,  the  user  can  then  vary  element  aspect  ratio,  the  arrangement 
of  elements  or  element  groups,  choose  uniform  or  variable  meshing,  and  even  introduce 
abrupt  mesh  changes.  All  of  these  will  influence  the  solution  to  some  degree. 

Constitutive  Model:  Assuming  an  appropriate  model  has  been  selected  to  account  for 
material  behavior  under  high-rate  loading,  criteria  for  material  failure,  and  descriptions  of 
post-failure  behavior,  the  problem  of  selecting  values  for  the  various  material  constants  in 
the  constitutive  model  remains.  These  must  be  selected  from  experiments  conducted  at  strain 
rates  appropriate  to  the  problem.  In  many  cases,  data  may  not  be  available.  This  is 
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particularly  true  for  situations  involving  material  failure.  The  estimation  of  these  factors  then 
depends  on  the  skill,  knowledge,  and  experience  of  the  user,  and  computed  results  will  vary 
accordingly. 

Contact  Surfaces  and  Material  Transport:  Lagrangian  codes  incorporate  a  wide  variety  of 
algorithms  to  account  for  contact-impact  situations.  Eulerian  codes  have  a  variety  of 
methods  for  determining  the  transport  of  material  from  one  cell  to  another.  Each  algorithm 
uniquely  affects  the  solution.  Some  codes  have  incorporated  a  variety  of  algorithms  and 
allow  the  user  a  choice.  The  burden  is  then  on  the  user  to  choose  wisely,  and  this  cannot  be 
done  without  a  knowledge  of  how  the  various  algorithms  affect  the  solution  of  both  global 
(e.g.,  displacements)  and  local  (e.g.,  strain)  variables. 

Shortcuts:  Because  of  the  expense  involved  in  3-D  calculations,  recourse  is  sometimes  made 
to  plane-strain  solutions.  These  are  2-D  approximations  of  3-D  phenomena.  Sometimes  they 
produce  excellent  results,  sometimes  disasters. 

This  report  focuses  on  problems  involving  meshing.  Subsequent  reports  will  address  the  other 
topics. 


5.  Problems  Related  to  Computational  Meshes 


Theoretically,  the  ideal  mesh  is  uniform  in  all  coordinate  directions  and  “converged”  for  the 
critical  variable  for  the  problem;  that  is,  it  is  small  enough  to  give  accurate  results  so  that  further 
refinement  dramatically  runs  up  the  cost  of  computing  with  negligible  improvement  in  accuracy. 
In  practice,  especially  when  performing  3-D  calculations,  this  goal  is  impossible  to  achieve. 
Compromises  must  be  made,  and  the  educated  analyst  must  know  what  effect  these  compromises 
have  on  the  numerical  solution  to  his/her  problem.  In  short  he/ she  must  know  the  physical  problem 
and  understand  how  compromises  to  the  computational  mesh  affect  the  solution. 
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There  are  a  number  of  factors  that  affect  mesh  integrity.  These  include  element  aspect  ratio, 
element  arrangement,  uniform  vs.  graded  meshes,  and  abrupt  changes  in  meshes.  Each  is  now 
looked  at  in  turn. 

5.1  Element  Aspect  Ratio.  Ideally,  calculations  would  be  done  with  1:1  aspect  ratios  in  the 
elements.  This  hardly  ever  happens  in  practical  2-D  and  3-D  calculations.  Thus,  it  would  be  nice 
to  know  how  the  solution  is  affected  when  elements  with  aspect  ratios  exceeding  1:1  are  used. 
Creighton  [20],  using  the  EPIC-2  code,  looked  at  the  effects  of  aspect  ratio,  artificial  viscosity,  and 
triangular  element  arrangement  on  elastic  and  elasto-plastic  impact  situations.  EPIC  uses 
constant-strain  triangular  elements  that  can  be  arranged  in  a  number  of  ways.  The  possibilities 
included  in  Creighton’s  study  are  shown  in  Figure  2.  The  elastic  calculations  were  performed  for 
a  steel  bar  with  length-to-diameter  (L/D)  ratio  of  100  striking  a  rigid  barrier  at  3.048  m/s.  Numerical 
solutions  were  compared  with  an  exact  solution  by  Skalak  [21],  Figure  3,  which  takes  into  account 
the  effects  of  radial  inertia.  Calculations  were  performed  with  one  (400  elements),  two 
(1,600  elements),  and  three  (3,600  elements)  crossed  triangles  (four  triangles  per  quadrilateral) 
across  the  radius  of  the  bar  (length,  L  -  12.7  cm;  diameter,  D  -  0.254  cm).  All  calculations  were 
done  with  an  aspect  ratio  of  1 : 1 .  The  axial  force  in  the  rod  as  a  function  of  position  was  compared 
with  Skalak’ s  solution  at  various  times. 

As  expected,  the  results  show  the  grid  acting  as  a  frequency  filter.  As  the  grid  is  resolved,  more 
and  more  ringing  is  computed  behind  the  wavefront.  Enough  high-frequency  components  are 
accounted  for  in  the  3,600-element  calculation  to  compare  very  favorably  with  the  analytical 
solution. 

Holding  the  number  of  elements  across  the  bar  radius  at  three,  the  element  aspect  ratio  was  now 
increased  from  1:1  to  4:1.  Again,  the  filtering  characteristics  of  the  mesh  is  clearly  shown. 
High-frequency  components  are  gradually  suppressed  until,  with  the  4: 1  mesh,  only  the  ringing 
directly  behind  the  wave  front  remains  (Figure  4). 
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(a)  JDJA=1  Hypotenuse  of  the  triangle  is  drawn  from  the 
lower  left  hand  corner  to  the  upper  right  hand 
corner. 


(b)  IDIA  =  2  Hypotenuse  of  the  triangle  is  drawn  from  the 
upper  left  hand  corner  to  the  lower  right  hand 
corner. 


(c)  Quadrilateral  elements  comprised  of  four  triar^les. 


Figure  2.  Element  Orientations  in  the  EPIC -2  Code. 
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Figure  3.  Comparison  of  Elementary  and  Skalak  Solutions  for  Bar  Impact  (Skalak  Solution 
Is  Drawn  Freehand). 
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Figure  4.  Signal  Filtration  With  Mesh  Aspect  Ratio. 


5.2  Element  Arrangement  The  arrangement  or  mixing  of  elements  within  a  computational 
mesh  can  be  a  problem.  Figure  2  shows  three  possible  arrangements  of  triangular  elements.  This 
aspect  of  mesh  generation  was  also  investigated  [20].  The  specific  problem  considered  was  the 
impact  of  a  steel  sphere  with  a  5.08-cm-thick  aluminum  target  at  a  velocity  of  1,524  m/s.  The 
two-triangle  orientations  (Figure  2[a]  and  [b])  gave  rise  to  asymmetries  in  the  calculation. 
Depending  on  the  orientation  of  the  diagonal,  results  were  either  too  stiff  or  too  soft.  Optimal 
performance  was  achieved  using  the  four-triangles-per-quadrilateral  grouping  (Figure  2[c]).  Similar 
results  had  been  obtained  by  Zukas  [22].  The  definitive  study  on  element  arrangement  for  accurate 
elasto-plastic  solutions  was  done  by  Nagtegaal  et  al.  [23]. 

Johnson  and  coworkers  [24, 25]  observed  that,  for  2-D  triangular  elements,  accuracy  is  much 
improved  if  the  elements  are  arranged  in  a  crossed-triangle  (four  triangles  within  a  quadrilateral) 
arrangement  or  if  a  single  average  pressure  is  used  for  each  group  of  two  adjacent  triangles  (two 
triangles  within  a  quadrilateral).  In  three  dimensions,  Johnson  et  al.  [25]  performed  a  number  of 
simulations  using  configurations  of  6  tetrahedral  elements  within  a  brick,  6  tetrahedral  elements  per 
brick  using  a  single  average  pressure  for  all  6  tetrahedra,  as  well  as  an  arrangement  of  24  tetrahedral 
elements  per  brick.  This  latter  case,  though  the  most  costly  configuration,  was  expected  to  produce 
the  best  results  since  the  asymmetries  present  in  the  six  tetrahedra  arrangement  would  be  absent.  It 
was  also  expected  that  the  six  tetrahedra  configuration  with  pressure  averaging  would  perform  better 
compared  to  the  six  tetrahedra  configuration  using  individual  element  pressures  due  to  the  reduction 
in  the  number  of  incompressibility  constraints  [23] .  Some  test  cases  were  compared  to  experimental 
data  [26].  Between  37,000  and  39,000  elements  were  used  for  these  calculations  of  an  L/D  -  10 
tungsten  rod  striking  a  2.5-cm  steel  plate  at  normal  incidence  with  a  velocity  of  1 ,520  m/s.  The  case 
of  the  six  tetrahedra  without  pressure  averaging  produced  the  closest  correlation  with  experiment  for 
residual  velocity  and  the  worst  for  residual  rod  length.  Furthermore,  there  was  indication  of  locking 
in  the  target  grid.  The  closest  correlation  for  rod  residual  velocity  and  residual  length  was  achieved 
for  the  symmetric  24  tetrahedra  per  brick  arrangement.  However,  in  all  cases,  including  a  2-D 
calculation,  EPIC  tended  to  overpredict  residual  length  by  about  16%.  This  could  be  due  to  a 
number  of  factors  not  necessarily  related  to  the  grid. 
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A  number  of  papers  present  comparisons  of  plastic  strain  profiles  for  Taylor  cylinder 
impacts — the  impact  of  a  deformable,  short  cylinder  (L/D  <  3,  generally,  although  longer  rods  have 
been  used)  striking  a  rigid  surface  at  impact  velocities  under  100  m/s.  Johnson  [24],  for  example, 
compares  DYNA,  NIKE,  and  EPIC-2  results  with  various  arrangements  of  triangles  (Figure  5). 
Many  more  papers  have  since  appeared  to  not  only  test  element  formulations  and  arrangements  but 
constitutive  relations  as  well.  Such  results  are  interesting  but  tend  to  show  only  minor  differences 
for  the  various  cases  considered.  It  is  now  clear  that  the  Taylor  cylinder  is  not  a  very  sensitive 
measure  of  element  arrangement  effectiveness.  Neither  is  it  a  good  discriminant  for  testing 
constitutive  models. 

5.3  Uniform  and  Variable  Meshes.  The  effects  of  mesh  size  on  wave  propagation  problems 
can  be  seen  in  Figures  6  and  7.  These  depict  the  impact  of  a  long  S-7  tool  steel  rod  (L/D  =10, 
hemispherical  nose,  D  =  1.02  cm)  into  2.56-cm-thick  RHA  plate  at  1,103  m/s.  The  characteristic 
lengths  of  this  problem  are  the  target  thickness  and  projectile  radius.  For  good  results,  there  should 
be  at  least  three  elements  across  the  radius  of  the  rod.  Since  uniform  mesh  spacing  is  the  ideal,  this 
governs  the  number  of  elements  to  be  used  in  the  projectile  and  target. 

Figure  6  shows  initial  grids  and  results  at  50  ms  after  impact.  One,  three,  and  five 
crossed-triangle  meshes  were  used  across  the  projectile  radius  for  the  coarse,  medium,  and  fine  cases, 
respectively.  The  target  grid  was  then  set  by  requiring  a  1:1  aspect  ratio  of  all  elements  in  the 
calculation.  The  coarse  grid  computation  shows  some  anomalies  near  the  projectile-target  interface 
and  a  v-shaped  crater,  indicative  of  numerical  difficulties  with  triangular  elements.  The  other  two 
calculations  (Figures  6[d]  and  [e])  appear  to  be  reasonable.  The  experimentally  determined 
projectile  residual  mass  was  32.1  g,  and  the  residual  velocity  was  690  m/s.  The  computed  residual 
masses  were  29  g,  36  g,  and  37  g  for  the  coarse,  medium,  and  fine  grids,  respectively,  while  the 
residual  velocities  were  600  m/s,  670  m/s,  and  680  m/s.  Coarse  zoning  is  adequate  for  “quick  and 
dirty”  scoping  calculations  or  to  get  preliminary  estimates  of  global  quantities  such  as  length  loss 
in  the  rod,  approximate  hole  size,  residual  velocity,  and  overall  deformed  shapes  of  the  two  solids. 
It  is,  however,  inadequate  to  resolve  strain  and  pressure  fields  with  any  degree  of  accuracy.  If  the 
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Figure  5.  Effects  of  Element  Arrangement  on  Plastic  Strain  Distribution  [24]. 
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Figure  7.  Gradual  vs.  Abrupt  Grid  Changes. 
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calculation  employs  failure  criteria  based  on  stress  wave  profiles,  results  from  such  coarse 
calculations  will  be  meaningless. 

As  spatial  resolution  is  increased  by  a  factor  of  2  for  2-D  simulations,  computer  time  rises  by 
a  factor  of  8  for  explicit  methods  [27].  Thus,  the  resolution  used  should  match  the  accuracy 
required.  If  comparison  is  to  be  made  with  time-resolved  pressure,  stress,  or  strain  data,  if  internal 
failure  (such  as  spall)  occurs,  fine  resolution  is  required  for  a  meaningful  calculation,  even  though 
the  cost  is  high.  Compromising  the  accuracy  of  the  calculation  to  save  money  in  these  cases  is  not 
justifiable  since  maximum  savings  come  from  not  doing  the  calculation  at  all.  This  has  the  added 
advantage  of  not  generating  meaningless  numbers  that  someone  not  familiar  with  the  events  of  the 
calculation  might  be  tempted  to  believe.  On  the  other  hand,  if  only  global  (integral)  data  are 
available  for  comparison  (residual  masses,  lengths,  hole  sizes,  and  deformed  shapes),  a  reasonably 
crude  and  inexpensive  calculation  can  be  done,  provided  its  interpretation  is  not  pushed  too  far. 

It  is  possible  to  take  advantage  of  the  localized  nature  of  impact  problems  when  setting  up 
computational  grids.  Typically,  the  high  strain  rates  and  the  high  pressures,  strains,  and  temperatures 
that  accompany  them  are  confined  to  a  narrow  region,  or  process  zone,  that  extends  about  three  to 
six  striker  diameters  from  the  impact  interface,  depending  on  striking  velocity.  Figure  7(a)  shows 
results  for  a  calculation  that  takes  advantage  of  this  information  by  localizing  fine  zones  within  three 
diameters  of  the  impact  interface  and  then  gradually  expands  element  size  in  regions  where,  at  most, 
elastic  waves  will  propagate.  The  result  is  accuracy  comparable  to  the  uniformly  gridded  case  but 
with  considerably  fewer  elements  and,  therefore,  considerable  savings  in  computer  time.  Care  must 
be  taken  not  to  change  element  sizes  too  rapidly,  however.  A  change  in  element  size  represents  a 
change  in  element  stiffness  even  if  material  properties  remain  the  same.  Traveling  waves,  when 
encountering  this  stiffness  difference,  will  act  as  if  an  impedance  mismatch  had  occurred.  Part  of 
the  wave  will  be  reflected,  part  transmitted.  If  element-to-element  size  variation  is  kept  below  1 0%, 
acceptable  results  can  generally  be  achieved.  Figure  7(b)  shows  the  negative  aspects  of  drastic 
changes  in  element  size. 
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Using  a  minimum  of  three  continuum  elements  across  a  critical  dimension  such  as  a  rod  radius 
or  plate  thickness  turns  out  to  be  a  good  rule  of  thumb  for  Lagrangian  calculations.  Eulerian 
calculations  require  considerably  more  [28] .  Some  problems,  however,  such  as  hypervelocity  impact 
or  self-forging  fragment  formation,  may  require  much  more.  These  are  situations  where  severe 
pressure  gradients  exist  and  move  with  time.  Also  situations  such  as  spall,  where  material  failure 
occurs  due  to  the  interaction  of  stress  waves  with  geometric  boundaries,  material  interfaces,  or  each 
other,  will  require  fine  resolution.  Zukas  et  al.  [29],  studying  explosively  formed  penetrator  (EFP) 
formation  and  penetration,  found  that  five  to  six  elements  through  the  liner  thickness  and  fine  zoning 
in  the  target  were  required  to  match  code  calculations  with  experiments.  Melosh  et  al.  [30]  modeled 
dynamic  fragmentation  on  a  laboratory  scale.  They  developed  a  fragmentation  model  and  found 
good  correlation  with  a  wealth  of  experimental  data  for  the  largest  fragment  size  and  the  fragment 
size-number  distribution,  provided  that  an  adequate  numerical  resolution  was  used.  Resolutions  of 
12  x  24  cells  (where  12  cells  defined  target  radii,  which  ranged  from  2-12  cm)  were  used  for  most 
calculations.  Resolutions  of  6  x  12  and  24  x  48  cells  were  also  used.  All  three  resolutions  gave 
about  the  same  result.  However,  substantially  finer  grids  (40  x  80)  were  needed  to  match  observed 
near-surface  spallation.  Johnson  et  al.  [25]  also  looked  at  the  effects  of  grid  size  on  fragment 
distribution  for  normal  impacts  of  copper  rods  at  2  km/s  against  1  -cm  steel  plates.  Calculations  were 
performed  using  1,600  (Case  A)  and  4,096  elements  (Case  B).  Figure  8  shows  the  effects  of 
gridding  on  fragment  size  distributions  for  the  different  grid  sizes. 

Johnson  and  Schonhardt  [31]  investigated  the  sensitivity  of  the  EPIC-2  and  EPIC-3  codes  to  a 
number  of  factors,  including  gridding  for  normal  and  oblique  incidence  problems.  Normal  impact 
calculations  for  an  L/D  -  10  tungsten  rod  striking  a  2.5-cm  steel  plate  at  a  velocity  of  1 .5  km/s  were 
made  with  1 44;  576;  1 ,296;  and  2,304  elements.  Oblique  impact  calculations  at  60°,  measured  from 
the  target  normal,  were  made  with  480;  1 ,920;  4,320;  and  7,680  elements.  With  the  exception  of  the 
lowest  resolution,  there  was  relatively  little  difference  for  the  higher  resolution  calculations  in  terms 
of  residual  velocity,  residual  mass,  and  hole  diameter  for  the  2-D  calculations.  Significant  increases 
in  CPU  times  were  observed  as  resolution  was  enhanced.  Similar  results  were  obtained  for  the  3-D 
simulations  (Figure  9)  comparing  residual  masses,  velocities,  hole  diameters,  rotational  velocity,  and 
deflection). 
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Figure  9.  Grid  and  Contact  Surface  Effects  for  Oblique  Penetration  [31]. 
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Zoning  requirements  for  2-D  and  3-D  Eulerian  calculations  with  the  CTH  code  for  long-rod 
penetration  problems  are  considered  in  Littlefield  and  Anderson  [28]. 

Guidelines  for  determining  resolution  and  material  modeling  for  practical  engineering  problems 
are  given  in  the  report  of  the  National  Materials  Advisory  Board  Committee  titled  “Materials 
Response  to  Ultra-High  Loading  Rates”  [16].  The  committee  recommended  an  iterative  procedure 
of  successive  refinements  involving  computations  with  existing  relatively  simple  failure  descriptions, 
dynamic  material  characterization  employing  relatively  simple  and  standardized  techniques,  and 
experimentation  to  produce  useful  results  for  design  purposes  in  many  applications.  Their  report 
suggests  that 

.  .  .  rough  computations,  using  simple  material  models  with  published  or  even 
estimated  material  properties,  may  be  used  in  conjunction  with  exploratory  test 
firings  to  scope  an  initial  design.  Comparison  of  test  data  with  predictions  may 
reveal  discrepancies  which  suggest  refinements  in  the  computations  or  material 
models,  and  the  need  for  some  material  property  measurements.  Once  reasonable 
agreement  has  been  achieved,  another  round  of  computations  may  then  be  performed 
to  refine  the  design.  Test  firings  of  this  design  might  use  more  detailed  diagnostic 
instrumentation.  This  sequence  is  iterated,  including  successively  more  detail  in 
computational  models,  material  property  tests,  and  ordnance  test  firings,  until  a 
satisfactory  design  is  achieved.  In  this  procedure,  unnecessarily  detailed 
computations,  material  property  studies  or  test  firings  are  minimized;  only  those 
details  necessary  to  achieve  a  satisfactory  design  are  included. 

5.4  Abrupt  Changes  in  Meshes.  Abrupt  grid  changes  are  common  in  statics  calculations.  The 
stress  gradient  is  fixed  in  space.  A  grid  sufficiently  fine  for  all  practical  purposes  is  superimposed 
over  the  region.  The  remainder  of  the  physical  body  is  then  modeled  with  rather  large  elements  to 
account  for  the  total  mass  and  boundary  conditions. 
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In  dynamics  calculations,  the  location  of  the  stress  or  pressure  gradient  is  a  function  of  time  as 
well  as  space.  It  cannot  be  overemphasized  that  wave  propagation  governs  the  response.  Recall 
from  wave  propagation  theory  that  stress  waves  are  reflected  from  material  interfaces,  geometric 
boundaries  and  each  other.  Partial  reflection  also  occurs  at  internal  points  of  discontinuity. 
Chapter  10,  in  Fried  [32],  describes  this  very  lucidly,  as  do  Bazant  and  his  colleagues  [33-35]  for 
a  number  of  practical  systems.  A  spurious  reflection 

. . .  occurs  when  the  traveling  wave  crosses  over  from  a  fine  mesh  region  to  a  coarse 
mesh  region.  A  given  mesh  size  has  a  lower  limit  to  the  wavelength  it  can 
approximate  or  an  upper  limit  to  the  frequency  it  can  transmit  -  the  cutoff  frequency. 

When  the  wave  enters  the  coarse  mesh  region,  some  of  its  high  frequency 
components  cannot  penetrate  and  are  reflected ....  As  the  wave  reaches  the  larger 
elements,  small  waves  appear,  traveling  backwards.  Also,  because  of  the  largest 
element  size,  dispersion  becomes  more  pronounced  in  the  form  of  leading  small 
waves  in  front  of  the  original  one  [32]. 

In  Fried  [32],  these  points  are  illustrated  with  examples  of  waves  in  strings.  Consider  also  the 
following  problem,  where  an  UD  ■=  3  mild-steel  projectile  impacts  an  armor-steel  target  at  a  velocity 
of  0.8  km/s  with  sharp  mesh  discontinuities  in  both  projectile  and  target.  Figure  10  shows  the  grid 
and  results  at  intermediate  times.  Figure  1 1  shows  the  end  result.  Note  that  the  grid  discontinuity 
has,  in  effect,  predetermined  the  outcome.  The  elements  are  sufficiently  large  that  the  distortions 
in  the  finely  gridded  zone  are  enhanced  by  reflected  high-frequency  waves  from  the  fine/coarse 
element  boundary.  Compare  the  results  of  Figure  1 2(a)  with  those  of  1 2(b),  where  a  graduated  mesh 
was  used,  and  Figure  12(c),  with  a  uniform  mesh.  Keep  in  mind  that,  in  high-velocity  impact 
problems,  the  most  severe  distortions  occur  within  three  to  six  projectile  diameters,  depending  on 
impact  velocity.  To  the  extent  possible,  this  region  should  be  uniformly  zoned,  with  element  size 
gradually  increased  by  no  more  than  10%  (folklore)  from  there  onward. 
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Figure  11.  Residual  Penetrator  and  Hole  Size  in  the  Presence  of  an  Abrupt  Grid  Change. 
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6.  Summary 


The  causes  of  disagreement  between  large-scale  code  calculations  and  reality  for  problems 
involving  fast,  transient  loading  have  been  discussed.  The  experience  of  the  analyst  is  a  prime  factor 
in  the  successful  use  of  commercial  codes  for  problems  in  dynamics.  The  analyst  needs  an 
appropriate  educational  background  and  a  keen  understanding  of  the  physics  and  mechanics  of  the 
problem  being  addressed.  He/she  must  also  have  sufficient  experience  with  numerical  techniques 
and  large-scale  computations  in  order  to  select  the  appropriate  computational  tool  for  the  problem 
and  to  evaluate  the  computational  results.  The  currently  available  commercial  codes  for  dynamics 
problems  are  in  no  way  “black  boxes”  that  can  be  assigned  to  junior  engineers  with  demands  for 
immediate  production.  The  consequences  of  inappropriate  analyses  can  range  from  embarrassment 
and  loss  of  funding  through  catastrophic  failure  of  poorly  designed  structures  under  service  loads, 
with  the  liabilities  and  litigations  that  inevitably  follow. 
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